function [] = runDynamicPrediction(runId)

rng(1); % for replicability

%%% Add paths
addpath('../general_funcs');
addpath('../model_funcs');
addpath('../model_funcs/dim_specific_funcs/dim2');
addpath('../model_funcs/dim_specific_funcs/dim3');
addpath('../consumerAdoptions');
addpath('../providerAdoptions');
addpath('../providerExits');

%%% Get input / output paths
config = getConfig_jointProcess(runId);

% Get config_consumerAdoptions and config_providerAdoptions: necessary input for make_Xdynamic_xxx functions
indivConfigs.consumerAdoptions = getConfig_consumerAdoptions(config.runId_consumerAdoptions);
indivConfigs.providerAdoptions  = getConfig_providerAdoptions(config.runId_providerAdoptions);
indivConfigs.providerExits      = getConfig_providerExits(config.runId_providerExits);

indivRunFolders.providerAdoptions  = config.runFolder_providerAdoptions;
indivRunFolders.providerExits      = config.runFolder_providerExits;
indivRunFolders.consumerAdoptions = config.runFolder_consumerAdoptions;

outputFolder = createFolderIfNotExist(sprintf('%s/dynamicPredictions', config.runFolder));
periodStartConditionsFolder = createFolderIfNotExist(sprintf('%s/periodStartConditions', outputFolder));

outputMatFile = sprintf('%s/Ypreds.mat', outputFolder);
outputMatFile2 = sprintf('%s/pred_vals.mat', outputFolder);

%%% Load data from the 3 submodels
data = load_data_fullmodel(indivRunFolders);

%%% Get initial conditions (from base case) + Xdynamic_staticInputs
[initialConditions, Xdynamic_staticInputs] = load_base_initial_conditions_and_Xdynamic_staticInputs(indivRunFolders);

%%%  Load parameter estimates from the 3 submodels
[~,paramsParts] = load_param_estimates(data, indivRunFolders);

%%% Aggregate the data from spell level to day level to get day-level covariates
data.consumerAdoptions = aggregateCtsTimeData_consumerAdoptions(data.consumerAdoptions);
data.providerAdoptions  = aggregateCtsTimeData_providerAdoptions(data.providerAdoptions);
data.providerExits      = aggregateCtsTimeData_providerExits(data.providerExits);

%%% Read dimensions
T  = data.consumerAdoptions.dims.dims1{4};
K1 = data.consumerAdoptions.dims.dims1{5};
K2 = data.consumerAdoptions.dims.dims1{6};
K  = data.providerAdoptions.dims.dims1{3};

%%% Combine dims into object mydims
mydims.T  = T;
mydims.K  = K;
mydims.K1 = K1;
mydims.K2 = K2;

% Get labels
periodLabels = data.periodLabels;
geoLabels    = data.geoLabels;

%%% Compute the parts of V that does not depend on installedBases
XbetaFixed = calc_Xbeta_fixed(data, paramsParts);

%%% Preprocess Xdynamic_staticInputs inputs to compute the parts of V that depend on installedBases
Xdynamic_staticInputs = permute_Xparts_static_X_X_FEs(Xdynamic_staticInputs, mydims);


%%%% Launch predictions in base-case scenario
NumPreds = 100;
Ypreds.consumerAdoptions = zeros(T,K,NumPreds);
Ypreds.providerAdoptions  = zeros(T,K,NumPreds);
Ypreds.providerExits      = zeros(T,K,NumPreds);
tic
for pp = 1:NumPreds
	[Ypred_pp, periodStartConditions_pp] = generate_fullmodel_sequence_ctsTime(indivConfigs, initialConditions, XbetaFixed, Xdynamic_staticInputs, paramsParts); % it is Y_tk1 that is going to be stored (too much storage required for Y_tk1tk2)
	Ypred_pp.consumerAdoptions = reshape(Ypred_pp.consumerAdoptions, [T*K1 K2]); % (T*K1) x K2
	Ypreds.consumerAdoptions(:,:,pp) = reshape(sum(Ypred_pp.consumerAdoptions, 2), [T K1]); % T x K1
	Ypreds.providerAdoptions(:,:,pp)  = Ypred_pp.providerAdoptions;
	Ypreds.providerExits(:,:,pp)      = Ypred_pp.providerExits;
	
	% Save random states in separate file for each pp
	periodStartConditions = periodStartConditions_pp;
	save(sprintf('%s/periodStartConditions_%d.mat', periodStartConditionsFolder, pp), 'periodStartConditions');
	
	% Display time
	displayRemainingTime(pp, NumPreds, 1);
end
toc

%%% Save predictions to output file 1
save(outputMatFile, 'Ypreds', 'periodLabels', 'geoLabels', '-v7.3');

% Get predictions of lagged base + aggregate predictions + save those to output file 2
pred_vals = compute_laggedBases(Ypreds, initialConditions, false);
pred_aggVals = calc_Yagg(pred_vals);
save(outputMatFile2, 'pred_vals', 'pred_aggVals', 'periodLabels', 'geoLabels', '-v7.3');


% Check consistency between pred_vals and periodStartConditions
clear periodStartConditions;
for pp = 1:NumPreds
	load(sprintf('%s/periodStartConditions_%d.mat', periodStartConditionsFolder, pp), 'periodStartConditions');
	for tt = 1:T
		if ~isequal(squeeze(pred_vals.laggedProviders(tt,:,pp))',              periodStartConditions{tt}.laggedProviders); error('1'); end;
		if ~isequal(squeeze(pred_vals.laggedConsumersOrig(tt,:,pp))',         periodStartConditions{tt}.laggedConsumersOrig); error('2'); end;
		if ~isequal(squeeze(pred_vals.popAtRisk.providerAdoptions(tt,:,pp))',  periodStartConditions{tt}.popAtRisk.providerAdoptions); error('3'); end;
		if ~isequal(squeeze(pred_vals.popAtRisk.consumerAdoptions(tt,:,pp))', periodStartConditions{tt}.popAtRisk.consumerAdoptions); error('4'); end;
	end
end

end
